#######################################################################
########REPLICATION FILE: RACIAL ISOLATION DRIVES RACIAL VOTING########
########################POLITICAL BEHAVIOR#############################
####################MELISSA SANDS, DANIEL DE KADT######################
#######################################################################

rm(list=ls(all=T))
library(foreign)

####################################################
######## ME PLOT FOR INDIVIDUAL SURVEY DATA ########
####################################################

##### Function for Handling Stata Output ##### 
numerical = function(x){ 
  as.numeric(
    unlist(strsplit(as.character(x), split='=', fixed=TRUE))[2]
  )
}

##### Set Working Dir & Call Data ##### 
setwd("YOUR DIRECTORY HERE")

d = read.dta("data\\rugplot_data.dta")

##### Subset Only To Those Who Voted ##### 
d1 = subset(d, d$vote==1)

##### Coefficients From Column 3 of Indvidual Analysis ##### 
modelresults = read.csv("results\\individual_maineffects.csv")

modelresults = apply(modelresults,  c(1, 2), numerical)

coef_maineffect <- modelresults[1,4] #white
coef_interactionterm <-  modelresults[7,4] #white_whiteiso
variance_maineffect <- modelresults[2,4]^2  #var white
variance_interactionterm <-  modelresults[8,4]^2 #var interaction
COVARIANCE_maineffect_interactionterm <- -.0012984 #cov white, white_whiteiso, calculated
degrees_of_freedom <-  1114  #dof from model

ME <- function(whitefrac){
  (coef_maineffect + coef_interactionterm * whitefrac)
}

SE <- function(whitefrac){
  sqrt(
    variance_maineffect +
      whitefrac^2*variance_interactionterm +
      2*whitefrac*COVARIANCE_maineffect_interactionterm
  )
}

df <-degrees_of_freedom 

confint.upper <- function(x){
  ME(x) - qt(.025,df)*SE(x)
}

confint.lower <- function(x){
  ME(x) + qt(.025,df)*SE(x)
}

pdf(file= "results\\white_white_iso_MEplot_ANCvs.pdf", width = 5, height = 5, pointsize = 10)  
plot(ME,  type="l" , ylim=c(-1,-.5), xlab="White Isolation 2011", ylab="Pr(ANC Vote)", axes=F)
par(new=TRUE)
plot(confint.upper, ylim=c(-1,-.5),  type="l", col="red",lty=2, xlab=NA, ylab=NA, axes=F)
par(new=TRUE)
plot(confint.lower,  ylim=c(-1,-.5), type="l", col="red",lty=2, xlab=NA, ylab=NA, axes=F)
legend(x=0.45, y=-0.6, bty="n", legend=c("Total Effect Estimate", "95% Confidence Interval", "Baseline White Effect"),
       lty=c(1,2,1), lwd=c(1,1,1), col=c("black","red","light grey")
       , y.intersp=0.8)
rug(d1$white_iso2011[d1$white==1],  side=1, col="blue")
rug(d1$white_iso2011[d1$white==0], side=3, col="red")
abline(h=coef_maineffect, col="light grey")
axis(1)
axis(2)
dev.off()
